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NOTATION 


Van  Driest  damping-length  parameter 
hydraulic  diameter 
dimensionless  stream  function 
u/u0 

reference  length  taken  equal  to  rQ 
local  Nusselt  number 
local  pressure 

pressure  at  inlet  to  entrance  section 

o 

dimensionless  pressure  drop  [=(p0  —  p)/(l/2)puQ] 
pressure  in  reservoir  upstream  of  inlet  plane 
Prandtl  number 
turbulent  Prandtl  number 

radial  distance  from  axis  of  revolution  {-rQ  -  y) 
pipe  radius  or  half  of  channel  width 
Reynolds  number  (=uQD/v,  uQL/v,  ugX/v) 

Reynolds  number  (=uee/v) 
local  Stanton  number 
temperature 

x-  and  y-component  of  velocity,  respectively 
centerline  velocity 

reference  velocity  taken  equal  to  averaged  velocity  in  duct 

shear  velocity  (= Aw/p) 

distance  along  body  surface 

dimensionless  axial  distance  (hxRqVd) 

distance  normal  to  x 
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centerline  distance 
boundary-layer  thickness 

OD 

displacement  thickness  [=^(1  - 
eddy  viscosity 

dimensionless  eddy  viscosity  (= 
transformed  y-coordlnate 
value  of  n  at  centerline 

OD 

momentum  thickness  [=^* (u/ue)(l 

parameter,  see  eq.  (4.5) 

intermittency 

Von  Karman's  constant 

kinematic  viscosity 

density 

shear  stress 
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l.o  INTRODUCTION 


The  purpose  of  the  present  paper  is  to  describe  an  efficient  finite- 
difference  method  for  the  solution  of  the  differential  equations  appropriate 
to  thin  shear  layers,  and  representing  conservation  of  mass,  momentum  and 
energy  and  to  demonstrate  that  calculations  of  internal  flows  and  flows  with 
separation  can  be  readily  obtained.  In  a  series  of  previous  papers,  see  for 
example  [1],  Cebeci  and  coauthors  have  reported  finite-difference  solutions 
of  boundary- layer  equations  for  a  range  of  boundary  conditions  corresponding 
to  external  boundary- layer  flows.  Essential  features  of  the  solution  method 
have  been  the  efficient  and  accurate  numerical  scheme  of  Keller  [2]  and,  for 
turbulent  flows,  a  simple  eddy-viscosity  formulation  [1],  Previous  applic¬ 
ations  have  been  mainly  concerned  with  steady  two-dimensional  and  axisymmetric 
flows  [1],  unsteady  two-dimensional  flows  [3],  steady  three-dimensional  flows 
[4]  as  well  as  with  the  solution  of  the  stability  equation  for  two-dimensional 
flows,  namely,  the  Orr-Sommerfeld  equation  [5],  More  recent  developments  have 
included  the  use  of  inverse  methods  which  allow  the  calculation  of  internal - 
flow  problems  [6,7,8]  and  calculations  through  small  regions  of  separated 
flow  [9,10,11].  As  a  consequence  of  this  previous  work,  the  greater  part  of 
this  report  is  concerned  with  the  method  described  in  reference  1  to  11  and 
Its  application  to  internal-  and  separated-flow  problems. 

A  small  number  of  calculated  laminar-flow  results  for  Internal  flows  are 
presented  and  discussed.  The  emphasis  is,  however,  on  turbulent  flows  and  on 
the  extent  of  present  ability  to  calculate  from  a  laminar  region,  through 
transition,  to  a  fully  turbulent  flow.  The  eddy-viscosity  formulation, 
previously  described  by  Cebeci  and  Smith  [1]  is  a  two-layer  algebraic  model. 

In  its  present  form,  it  cannot  be  expected  to  represent  three-dimensional 
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duct  flows,  but  it  is  shown  to  allow  the  accurate  calculation  of  developing 
turbulent  two-dimensional  plane-  and  round-duct  flows.  The  energy  equation 
is  solved  by  the  same  numerical  scheme,  and,  with  the  addition  of  a  turbulent 
Prandtl  number,  allows  correspondingly  accurate  results  for  heat  transfer 
properties.  Comparisons  with  measurements  of  heat-transfer  information  in  a 
duct  flow  which  undergoes  transition  from  laminar  to  turbulent  flow,  confirms 
the  validity  of  the  model  for  this  purpose. 

Previous  calculations  of  two-dimensional  duct  flows  have,  or  course,  been 
reported.  Developing  laminar  flows  have  received  considerable,  and  undeserved 
attention,  as  is  indicated  in  section  IV  of  this  report  and  corresponding  solu¬ 
tions  to  the  energy  equation  have  been  reported,  see  for  example  Tien  and 
Pawelek  [12].  Developed  turbulent-duct  flow  has  also  been  considered  in  detail, 
for  example,  Na  and  Habib  [13,14],  and  has  successfully  made  use  of  a  form  of 
the  present  turbulence  model.  More  recently,  attention  has  turned  to  the  solu¬ 
tion  of  the  momentum  and  energy  equations  for  developing  turbulent  duct  flow, 
for  example.  Nelson  and  Pletcher  [15],  Bradshaw,  Dean  and  McEligot  [16],  and 
Emery  and  Gessner  [17].  The  present  results  demonstrate  that  two-dimensional 
duct  flows  in  general,  and  including  those  which  pass  through  transition  from 
laminar  to  turbulent  flow,  can  be  represented  by  the  present  calculation  pro¬ 
cedure  which  combines  accuracy  with  efficiency.  Results  are  also  presented  for 
laminar  and  turbulent  separated  flows.  It  is  well  known  that,  for  a  prescribed 
external  velocity  distribution,  the  boundary-layer  equations  are  singular  at 
separation:  they  are  not  singular,  however,  when  the  displacement  thickness 
is  prescribed.  This  was  demonstrated  by  Catherall  and  Mangier  [18]  who  solved 
the  laminar  boundary-layer  equations  in  the  usual  way  until  separation  was 
approached  and  then,  by  assuming  a  displacement-thickness  distribution. 


calculated  the  external -velocity  distribution  and  local  flow  properties 
through  the  recirculation  region.  More  recent  investigations  of  the  inverse 
boundary-layer  method  have  been  reported,  for  example,  by  Klineberg  and 
Steger  [19],  Horton  [20],  Carter  [21,22],  Carter  and  Wornom  [23]  and  Williams 
[24,25],  and  have  involved  solutions  to  the  laminar  boundary-layer  equations 
with  prescribed  displacement-thickness  and  shear-stress  distributions  in 
regions  of  negative  shear  stress.  The  present  results,  also  obtained  with 
the  numerical  procedure  of  references  1  to  10,  are  for  laminar  and  turbulent 
flows  where  the  external  velocity  distribution  is  known.  They  correspond  to 
the  laminar  separating  and  reattaching  flows  previously  calculated  with 
Navier-Stokes  equations  by  Briley  [26]  and  Carter  [22]  and  to  the  near- 
separating  turbulent  boundary-layer  flows  of  Schubauer  and  Spangenberg  (see 
Coles  and  Hirst,  ref.  27).  The  results  again  demonstrate  that  the  procedure 
of  references  1  to  14,  modified  to  operate  in  an  inverse  manner,  allows 
accurate  results  with  very  small  computer  times. 

The  report  has  been  prepared  in  three  main  sections.  The  following 
section  states  the  basic  equations  and  boundary  conditions  corresponding  to 
the  present  range  of  flows;  the  equations  are  presented  first  in  physical 
coordinates  and  then  in  transformed  coordinates,  appropriate  to  boundary-layer 
flows.  Numerical  methods  are  reviewed  in  the  first  subsection  of  Section  3 
and  alternative  approaches  to  internal -flow  problems,  making  use  of  the 
numerical  scheme  of  references  1  to  14,  are  described  in  subsections  3.2 
and  3.3.  The  nonlinear  eigenvalue  approach  of  section  3.2  is  simpler  than 
the  alternative  Mechul  function  approach  of  section  3.3  and  is  preferred 
for  internal-flow  calculations.  The  Mechul  method  is,  however,  also  appropriate 
to  internal  and  external  flows  with  separation,  and  its  application  to 
separated  external  flows  is  described  in  subsection  4.5.  The  solution  of 
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the  energy  equation,  for  internal  flows,  is  considered  in  section  3.4. 

The  results  are  presented  in  section  4  which  emphasizes  the  capabilities 
of  the  eddy-viscosity  formulation  for  turbulent  flows  and  of  the  numerical 
schemes.  A  closing  section  presents  brief  conclusions. 


I 


2.0  BASIC  EQUATIONS  AND  BOUNDARY  CONDITIONS 


The  equations  used  to  characterize  the  flow  configurations  of  section  4 
may  be  written  in  the  form 


Continuity 


Momentum 


37  +  |y  <rk,l  ■  0 


(2.1) 


Energy 


„  3U  .  „  3U  .  I  d£.  v  3  r^kn  4.  9ul  lo  o\ 

u  a7+  v  37'"^d7+7^[r  0  +  V  ayj  (2-2) 

3T  .  aT  _  v  1  a  f  k  /.  Pr  +\  aT  T  ,9  ,, 

u  a7+ v  a7‘  PF^iy  |.r  l1  +  P^em)ayJ  (2-3) 


They  represent,  respectively,  conservation  of  mass,  momentum  and  energy  for 
two-dimensional  boundary-layer  flows  where  variations  in  density  and  specific 
heat  are  unimportant.  For  turbulent  flows,  they  presume  that  the  diffusion  of 
momentum  and  energy  may  be  represented  by  an  effective  viscosity  and  Prandtl 
number,  respectively:  the  specific  assumptions  used  here  are  introduced  and 
discussed  in  section  4.  The  flow  index,  k,  is  zero  for  plane  flow  and  unity 
for  axi symmetric  flow. 

The  parabolic  nature  of  the  thin  shear-layer  (TSL)  equations  requires  that 
boundary  conditions  be  provided  on  three  sides  of  the  solution  domain.  In 
addition  to  initial  conditions  for  the  momentum  and  energy  equations,  the  zero 
slip  condition  requires  that  the  u-velocity  be  zero  at  the  wall  and  the  wall 
temperature  or  heat-flux  may  be  prescribed  according  to  the  problem  to  be 
solved. 
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Equations  (2.1)  to  (2.3)  may  be  solved  in  the  forms  presented  above  or 
may  be  transformed,  for  a  particular  problem,  to  allow  solutions  with  improved 
economy  and  accuracy. 

2.1  Internal  Flows 

For  Internal  flows,  where  a  symmetry  condition  usually  provides  the  third 
boundary  condition,  it  is  convenient  to  transform  the  equations  to  the  Falkner- 
Skan  form  [28]  to  obtain  solutions  in  regions  where  the  boundary-layer  thick¬ 
ness  Is  less  than  half  the  duct  width.  Where  the  boundary-layer  fills  the 
duct,  physical  coordinates  are  desirable.  In  the  present  calculations,  the 
change  from  transformed  to  physical  coordinates  is  made  when  the  boundary-layer 
thickness  becomes  approximately  0.75  of  the  half-width  of  the  duct. 

It  is  convenient  to  nondimensionalize  equations  (2.1)  to  (2.3)  and,  for 
solutions  in  transformed  variables,  to  introduce  Mangler's  transformation  and 
the  stream  function.  The  introduction  of  the  dimensionless  quantities 

p*  =  “T  »  u*  =  jp  »  v*  *  J-  y*  -  £  x*  =  £  ,  r*  *  f- 
puQ  uo  o  L  L  L  L 

uL  j  (2-4) 

RL  "  ~  *  9  "  T" 

o 

allows  equations  (2.1)  to  (2.3)  to  be  expressed  as: 


|jr[(r*)ku*]  +  fyrt(r*)kv*]  .  0 

(2.5) 

U*  3U*  +  *  3U*  .  42^  +  1  8  (b 

ax*  v  ay*  dx*  (r*jk  ay*  \bl 

au*\ 

ay  *1 

(2.6) 

u*  +  v*  lfl_  =  J _ L _L_/b  1SL\ 

3x*  v  ay7  Pr  (r*jk  ay*  \b2  ay*) 

(2.7) 

6 


where 


b,  «  (r*)k(l  +  e*),  b2  =  (r*)k(l  +  Pr/Prt  e*)  (2.8) 

The  Mangier  transformation 

dx  -  (r*)2kdx*,  dy  =  (r*)kdy*  (2.9) 

ensures  that  equations  (2.5)  to  (2.7)  have  a  form  which  is  common  to  plane 
and  axi symmetric  flows,  i.e. 
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Here  u  and  v  are  related  to  u*  and  v*  through 


U*  =  u,  v*  - - 2_j_  v - ^-r-  4^ 

(r*)k  (r*)k  d> 

Introducing  the  stream  function,  defined  as 


u  dy 


5-4. 


V  =  — 


equations  (2.10)  to  (2.12)  may  be  rewritten  as 


(b  F")'  =  421  +  f>  ILL  -  F"  — 

1  dx  3x  3X 

(b~g ' ) *  -  Pr  (f1  4-9’ 

C  \  3X  3 X  / 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


Here  primes  denote  differentiation  with  respect  to  y. 
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Equations  (2.15)  and  (2.16)  may  be  rewritten  in  the  transformed  variables 
x  and  n  by  introducing  the  Falkner-Skan  variable. 


n  =  yfiji 


and  the  modified  stream  function 


Thus 


F(x,y)  =  f(x.n) 


(b,f") '  +  l  ff"  =  x  (&-  +  f'  _  f 

\dx  3x  3X  / 

(b2g')'  +  1  Pr  fg*  =  Pr  x  If'  ^.-g'  ^ \ 

\  3X  3 X  / 


(2.17a) 


(2.17b) 

(2.18) 

(2.19) 


where  the  primes  now  represent  differentiation  with  respect  to  n.  The 
parameters  b.|  and  b2  may  also  be  rewritten  with  r*  given  by 


(2.20) 


When  the  governing  equations  are  solved  in  transformed  variables,  the 
boundary  conditions  are  specified  at  the  wall  and  at  the  edge  of  the  boundary 
layer.  For  the  case  of  no  mass  transfer,  they  are  given  by: 


y  3  0,  u  =  0,  v  ®  0,  aT  (x)  +  bT'(x)  =  given 

rf  W 


y  =  6, 


u  =  ue. 


T  *  T 


(2.21a) 

(2.21b) 


Here  a  =  0  corresponds  to  the  case  in  which  wall  heat  flux  is  specified  and 
b  =  0  corresponds  to  specified  wall  temperature.  When  the  governing  equations 


I 
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are  solved  in  physical  variables,  the  boundary  conditions  are  specified  at  one 
wall  and  at  the  centerline  of  the  duct  provided  both  surfaces  are  kept  at  the 
same  temperature.  Otherwise  the  centerline  boundary  conditions  must  be 
replaced  by  those  on  the  other  wall,  similar  to  those  in  (2.21a)  except  with 
different  Tw(x).  For  the  case  when  we  use  the  boundary  conditions  on  the 
centerline  y  =  yc, 

y  =  V  ly’0'  f^*°  <2-22> 

In  terms  of  dimensionless  quantities,  equations  (2.21)  and  (2.22)  can  be 
written  as 


Wall  boundary  conditions 


n  =  0 

f  =  f‘  =  0, 

+  bg^(x)  =  given 

(2.23a) 

y  =  o 

O 

II 

II 

u. 

+  bg*(x)  =  given 

(2.23b) 

Boundary-layer  edge  conditions 


Centerline  boundary  conditions 

y  =  ^  f"  *  0  g‘  =  0  (2.25) 

The  presence  of  the  dp*/dx  term  in  either  equation  (2,15)  or  (2.18) 
introduces  an  additional  unknown  to  the  system  given  by  (2.15),  (2.16),  (2.23b) 
and  (2.24)  or  by  (2.18),  (2.19),  (2.23a)  and  (2.25).  Thus,  another  equation 
is  needed  and  is  provided  by  the  conservation  of  mass  In  Integral  form.  In 
the  case  of  a  two-dimensional  duct  with  constant  cross  section,  mass  balance 
gives 


Here  L  corresponds  to  the  half-width  of  the  duct.  In  terms  of  dimensionless 
variables  and  stream  function  expressed  in  physical  coordinates,  equation 
(2.26)  can  be  written  as 


F(x,^p  =  ^  (2.27) 

When  the  stream  function  is  expressed  in  transformed  variables,  equation  (2.26) 
can  be  written  in  a  form  similar  to  equation  (2.27),  i.e. 

f(x,nL)  =  <^[/x  (2.28) 

Similarly,  for  the  case  of  a  circular  pipe,  relations  similar  to  (2.27)  and 
(2.28)  can  be  written  as 

F(x,^)=^v1^  (2.29) 

and 

f(x,nL)  =  j  *^[7x  (2.30) 

2.2  Separated  Flows 

As  is  well  known  ,for  a  given  external  velocity  distribution 
the  governing  TSL  equations  become  singular  at  separation.  This  occurs  when 
wall  shear  is  equal  to  zero.  To  continue  the  calculations  into  the  region  of 
negative  wall  shear,  It  is  necessary  to  prescribe  an  alternative  boundary 
condition,  such  as  wall  shear  or  displacement  thickness  and  to  treat  the 
external  velocity  or  pressure  as  an  unknown.  According  to  the  Mechul  method, 
to  be  described  In  section  3.3,  separated  flow  calculations  for  external  flows 
can  be  achieved  by  solving  the  TSL  equations  subject  to  the  following  boundary 
conditions: 
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3.0  NUMERICAL  METHODS 


The  purpose  of  this  section  is  to  describe  numerical  methods  for  the 
solution  of  the  parabolic,  thin-shear-layer  form  of  the  momentum  and  energy 
equations.  In  the  following  subsection,  a  brief  review  of  appropriate 
numerical  methods  is  presented  and  the  relative  advantages  of  implicit  and 
explicit  finite-difference  methods  is  discussed.  As  a  result,  the  Box  scheme 
Is  preferred  and  is  described  in  detail  in  relation  to  the  momentum  equation. 
The  application  of  the  Box  scheme  to  the  solution  of  the  internal  flow  problem 
of  interest  here  is  described  in  the  following  two  subsections  which  are  con¬ 
cerned  with  a  "nonlinear  eigenvalue"  and  the  "Mechul"  approach,  respectively: 
both  methods  were  developed  by  Keller  and  Cebeci,  see,  for  example  references  6 
and  9.  In  the  former  approach,  the  pressure  is  regarded  as  an  eigenvalue 
parameter  and  in  the  latter  as  an  unknown  function.  As  indicated  in  reference 
9,  the  Mechul  approach  has  the  advantage  of  allowing  solutions  in  regions  of 
flow  recirculation  but  the  nonlinear  eigenvalue  method  is  simpler  to  use.  The 
penultimate  subsection  describes  the  solution  of  the  energy  equation  and  the 
section  ends  with  brief  corrments  on  the  solution  algorithm. 

3.1  A  Brief  Review  of  Numerical  Methods  for  TSL  Equations 

Various  methods  for  solving  TSL  equations  have  been  described  in  reference 
28  and  Include  finite-difference  methods,  method  of  lines,  Galerkin  method  and 
finite-element  methods.  Of  these  methods,  the  finite-difference  methods  are 
much  more  appropriate  to  solve  the  TSL  equations  and,  as  a  consequence,  only 
finite-difference  methods  are  considered  here. 

To  illustrate  the  application  of  finite-difference  methods  to  parabolic 
equations,  it  is  useful  to  consider  the  unsteady,  one-dimensional  equation, 
appropriate  to  the  Rayleigh  problem  and  to  heat  conduction,  i.e. 
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/ 


3U  _  3  U 

9t '  17 


0  <  x  <  L 


(3.1) 


The  solution  of  (3.1),  as  in  any  partial -differential  equation,  requires 
initial  and  boundary  conditions,  which  serve  to  complete  the  formulation 
of  a  "meaningful"  problem.  Let  us  assume  that  they  are  given  by 


Boundary  conditions 


Initial  conditions 


x  =  0 
x  =  L 


t  =  0 


u  =  given 
u  =  given 


u  =  u(x)  =  given 


(3.2a) 


(3.2b) 


The  solution  of  (3.1)  by  a  finite-difference  method  requires  that  we 
approximate  the  derivatives  by  some  difference  approximations.  Using  Taylor's 
theorem  and  requiring  that  the  function  u  and  its  derivatives  are  single¬ 


valued,  finite  and  continuous  functions  of  x,  we  can  write 

u(x  +  h)  =  u(x)  +  hu'(x)  +  j  h2u"(x)  +  g-  h3u'"  (x)  + 
and 

u(x  -  h)  =  u(x)  -  hu' (x)  +  h2u"(x)  -  g-  h3u"' (x)  + 

Adding  (3.3a)  and  (3.3b)  we  get 

u(x  +  h)  +  u(x  -  h)  =  2u(x)  +  h2u"(x) 
provided  the  fourth  and  higher-order  terms  are  neglected.  Thus, 

u"(x)  =  ^5-  [u(x  +  h)  -  2u(x)  +  u(x  -  h)] 
hz 
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with  an  error  of  order  h  . 

If  (3.3b)  is  subtracted  from  (3.3a) 

u'(x)  •  Irr  [u(x  +  h)  -  u(x  -  h)] 


(3.3a) 


(3.3b) 


(3.4) 


(3.5) 


Equation  (3.5)  approximates  the  slope  of  the  tangent  at  P  by  the  slope 
of  the  chord  AB  (see  figure  1)  and  is  called  a  central -difference  approxima¬ 
tion.  The  slope  of  the  tangent  at  P  can  also  be  approximated  by  either  the 
slope  of  the  chord  PB,  giving  the  forward -difference  formula 

u'(x)  =  jj-  [u(x  +  h)  -  u(x)]  (3.6) 

or  the  slope  of  the  chord  AP  giving  the  backward-difference  formula 

u'  (x)  =  j-j-  [u(x)  -  u(x  -  h)]  (3.7) 

Note  that  while  (3.5)  is  of  0(h2),  both  (3.6)  and  (3.7)  are  of  0(h). 

To  demonstrate  finite-difference  notation,  figure  2  indicates  a  set  of 
net  points  on  the  x-t  plane,  i.e. 

xi  =  xo  +  i h *  i  =  0,  1,2,  ...  I 

and  (3.8) 

t.  =  tQ  +  jk,  j  =  0,  1 ,  2,  ...  J 


Figure  1.  Notation  for  approximation  Figure  2.  General  finite-difference 
of  derivatives.  grid  notation. 


where  h  and  k  are  the  distance  and  time  mesh  width,  respectively.  With 


u(tj,xi)  =  (3.9) 

the  difference  approximations  of  equations  (3.4)  to  (3.7)  may  be  written  as 


u'(x^) 


1 

2F 


[uJ  -  — 


i+1 


ui-i] 


*  < 


i+1  ui] 


Ji-1 


] 


central  difference  from  (3.5) 
forward  difference  from  (3.6)  (3.10) 

backward  difference  from  (3.7) 


u"(x^) 
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H.i  -2ui  *  ui-i] 


(3.11) 


The  solution  of  (3.1)  may  be  obtained  by  using  either  an  explicit  or  an 
implicit  method.  In  the  former,  the  value  of  u  at  the  downstream  station  is 
expressed  in  terms  of  upstream  quantities  and  the  corresponding  equation  can  be 
solved  explicitly:  in  the  latter,  u  at  the  downstream  station  is  expressed 
in  terms  of  its  downstream  neighbors  and  the  known  upstream  quantities. 

An  explicit  formulation  may  be  obtained  by  representing  su/9t  by  the 
forward-difference  example  of  (3.10)  and  3^u/ax^  by  (3.11),  both  at  the  net 
point  (x-j,  tj)  as  indicated  in  figure  3;  i.e. 

UJ+1  -  u^ 

~  -k--1  •  rr  Cui*i  - 2ui  *  “i-i 3  <3-12> 

or  h 

ui+1  =  ui  +  [ui-l  “2ui  +  ui+l]’  1-1.2,  3,  ...I  (3.13) 

Equation  (3.12)  provides  an  explicit  finite-difference  approximation  to  (3.1) 
centered  at  the  net  point  (x-j,  tj).  The  value  of  u^+1  is  expressed  in  terms  of 


o  UNKNOWN 
x  KNOWN 
O  "CENTERING" 


Figure  3.  Finite-difference  grid  for  an  explicit  method. 

upstream  values  which  are  known  and  the  equation  allows  the  value  of  u  to 
be  obtained  at  all  values  of  tj+-|.  The  numerical  error  inherent  in  this 
scheme  can  be  shown  to  be  of  order  (k  +  h2)  and,  as  a  result  the  x-step,  h, 
must  be  maintained  small  to  ensure  acceptable  accuracy.  In  addition,  and 
although  explicit  formulations  are  computationally  simple,  they  can  lead  to 

instabilities  unless  the  time  step  is  also  small;  in  this  connection  it  is 
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necessary  to  ensure  that  k  <  1/2  h  . 

In  contrast,  implicit  methods  do  not  have  stability  requirements  and 
significantly  larger  time  steps,  with  corresponding  economy,  can  be  taken. 

As  an  example,  the  Crank-Nicholson  [29]  approach  replaces  3^u/8x  by  the 
means  of  the  finite-difference  approximations  at  (x..,  tj_y2);  i.e. 


and,  with  the  central -difference  form  of  (3.10)  (figure  4)  leads  to 


Figure  4. 


□  UNKNOWN 
x  KNOWN 
O  "CENTERING" 


t 


J 

j+1 1 

j5 

1 

! 

1— 

t— 

< 

) 

| _ 

[■■■  !  f 

1 - 1 

i  i+1 

Finite-difference  grid  for  Crank-Nicholson  method. 


u^+1  -  J 

1  i 


1 

2 


—  2u^+^  +  u^]  u*  ,  —  2U1?  +  u^  , 

1  l-l  .  i+l  l  l-l 

-  -  +  - -jj - 


Equation  (3.14)  can  be  rewritten  in  the  form: 


aiui-l  *  Vf’  *  ciuFl  ‘  r1 


where 


1  <  i  <  1-1 


ai 


7 


7 


b< =  *(’ +  7} 


Ci 


k 

7 


(3.14) 


(3.15) 


2  2 

The  errors  in  this  carefully  centered  scheme  are  of  order  (h  +  k  )  but  k 
need  not  be  related  to  h  for  stability  purposes.  The  scheme  is  uncondition¬ 
ally  stable  and  second-order  accuracy  may  be  achieved  with  uniform  x-spacing. 

On  the  other  hand,  the  unknown  value  of  u  is  expressed  in  terms  of  five  other 
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values  of  u  with  two  of  them  unknown:  as  a  consequence,  the  computational 
arithmetic  is  more  extensive. 

The  solution  of  (3.15)  can  be  obtained  very  easily.  In  vector  notation, 
the  equation  may  be  written  as 

A  Jd  =  J:  (3.16) 


where  we  have  introduced  the  (I  —  1 )  dimensional  vectors 


"U1 

"ri  " 

u2 

r2 

r  = 

• 

• 

% 

UI-1 

m  M 

rI-l 

_ 

(3.17a) 


and  the  (I-l)-order  matrix  (called  tridiaqonal  matrix) 


(3.17b) 


Then  the  solution  of  (3.16)  is  obtained  by  two  sweeps.  In  the  so-called 
forward  sweep,  the  followina  are  determined; 


61 

=  v 

Yl  -  C,/#, 

6i 

=  bi  "Vi-1 

1  *  2,  3 . 1-1 

Yi 

=  C1/B1 

1  =  2,  3,  ....  1-2 

(3.18a) 
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In  the  backward-sweep. 


1 


Z1  =  rl/8l*  Z1  =  (ri  -aizi-l)/Bi  1  -  2,  3 . 1-1 

(3.18b) 

UI-1  =  ZI-1  ui  =  zi  ~Wl  i  =  1-2,  1-3 . 1 

are  obtained. 

An  alternative  implicit  method  due  to  H.B.  Keller  [2]  and  Keller  and 
Cebeci  [ 6 J  is  now  described  and  is  referred  to  as  the  Box  method.  This  method 
has  several  very  desirable  features  that  make  it  appropriate  for  the  solution 
of  all  parabolic  partial-differential  equations  as  discussed  by  Cebeci  and 
Bradshaw  [28]  and  Keller  [30].  The  main  features  of  this  method  are: 

1.  sliohtly  more  arithmetic  to  solve, 

2.  second-order  accuracy  with  arbitrary  (nonuniform)  x-spacing, 

3.  allows  very  rapid  x-variations. 

The  solution  of  a  parabolic  partial-differential  equations  by  the  Box 
method  require  the  followinq  four  steps: 

1.  reduce  the  aovernina  equations  to  a  first-order  system, 

2.  write  difference  equations  using  central  differences, 

3.  linearize  the  resulting  alqebraic  equations  (if  they  are  nonlinear) 
and  write  them  in  matrix-vector  form, 

4.  solve  the  linear  system  by  the  block  tridiagonal  elimination  method. 
Reverting  to  the  notation  cf  section  2,  equations  (3.1)  and  (3.2a)  may  be 

written  as: 


au  _  a  u 

,x '  ^7 


0, 


u  =  a; 


L, 


(3.19) 

(3.20) 
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-*■'  -  - A 


(3.21a) 


To  allow  (3.19)  to  be  written  as  a  first-order  system, 

u'  =  v 


v.  =  lit 


(3.21b) 


and  the  difference  equations  may  be  written  to  approximate  (3.20)  with  (3.21a) 


centered  at  P-^  (xj,yj_i/2)  (see  the  rectangle  of  figure  5)  to  aet; 


u  -  u 


1  =  Vj-l/2 


(3.22) 


a  UNKNOWN 
x  KNOWN 
O  "CENTERING" 


i-1  i 


Figure  5.  Finite-difference  grid  for  the  Box  scheme. 

The  difference  form  of  (3.21b)  can  be  obtained  by  centering  at  the  middle  of 
P1P2P3P4  (xi -l /2 »  yj-1/2)  as  Allows.  Center  (3.21b)  at  x i _ i /2  to  obtain 

IVllMV)1'1  ,  j.W-1  (3.23a) 


Equation  (3.23a)  can  also  be  written  as 


(v.)1  2  ui=.(v.)i-l  2  ui-l  iR1-l 


(3.23b) 
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In  the  second  step  equation  (3.23b)  is  centered  at  yj_i/2  to  obtain 

Here 


(3.23c) 


R1'!  =  . 

J-l/2 

•  ,  u1’1  +  uM 

u’-1  =  J _ Jzl 

Uj-l/2  2 


Vi  1  vi-l  2  i-1 

■V"  ^  j-i/2 


ill 


j -1/2 


Since  (3.23c)  is  linear,  as  are  the  corresponding  boundary  conditions,  the 
system  may  be  written  in  matrix-vector  form  as  shown  below. 


uo 

vo 

Uj 

Vj 

u 

\ 

b.c.  - - 

r  v " 

’  o  | 

0 

i 

i 

-nn 1 

■h 

i 

eq. (3.24c) 

•  -i 

i 

O  I 
2  • 

1 

0 

~T 

1 

| 

eq. (3.23) 

i  k. 

i 

- 

i 

-1  i 

k. 

1 

j_ . 

i 

r 

i 

0 

1  1 
| 

1 

i 

i 

1 

1 

1  0 

i 

0  i 

-1 

£ 

1 

i 

l 

1 - 

_ 1 

_  _ 

L. 

0  0 


0  0 


1  0 


r(riV 

((r2>W 


'(r2>y 


f(ri  V 

^r2*0/ 


(3.24) 


where 


‘■Vo’* 

<r2  )j  "  0 


(lVj  =  hj-lRj-l/2 

0  <  j  <  J-l 


1  iJiJ 


(r2)0  «  b 


This  equation  may  be  rewritten  as: 


*  im  l 


(3.25) 


where 


An 

0  0 

V 

\o~ 

B1  A1  C1 

\  \  \ 

$1 

• 

• 

• 

JCi 

Bj  Aj  cj 

\  \  \ 

BJ-1  AJ-1  CJ-1 
\  \a 

6  = 

'V 

JC  * 

BJ  AJ 

ij 

(3.26) 


(3.27a) 


and  A.,  B  ,  C.  are  2x2  matrices  defined  as  follows: 

J  J  J 


Note  that,  as  in  the  Crank  Nicholson  method,  the  implicit  nature  of  the 
method  has  again  generated  a  tridiagonal  matrix,  but  with  2x2  blocks  rather 
than  scalars. 
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The  solution  of  (3.25)  is  obtained  by  the  block  elimination  method  dis¬ 
cussed  by  Keller  [31]  and  by  Cebeci  and  Bradshaw  [28].  This  method  consists 

of  two  sweeps.  In  the  forward  sweep  we  compute  r.-,  A;  and  w.  from  the 

J  J 

following  recursion  formulas 


“o  =  Ao 

Yj-i =  Bj 

AJ  ■  Aj  -  rjcj-l 

«o  =  to 

■  tj  -  rJ«j-l 


1  iJiJ 


1  i  j  <  J 


Here  has  the  same  structure  as  B^,  that  is. 


'(“n  (ai2>A 

V  0  0  / 


and  the  second  row  of  a-  is  the  same  as  the  second  row  of  A., 

J  J 

/'“ll’j  (o12>A 


(3.28) 


In  the  backward  sweep,  is  computed  from  the  following  recursion  formulas: 


4j4) =  ( 

“jij  *  ~  J  '  J-'.  J-2 . o 

For  further  details,  the  reader  is  referred  to  Cebeci  and  Bradshaw  [28]. 


(3.29) 


3.2  Momentum  Equation  and  the  Nonlinear  Eigenvalue  Method  for  Internal  Flows 
In  the  mathematical  description  of  an  internal  boundary-layer  flow,  the 
geometry  rather  than  the  pressure  distribution  is  known.  Thus,  the  solution 
of  (2.18)  requires  a  form  of  iteration  to  determine  dp*/dx.  In  the  nonlinear 
eigenvalue  approach,  this  situation  is  performed  by  applying  Newton's  method 
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to  the  solution  of  the  nonlinear  algebraic  equations  corresponding  to  a 
finite-difference  form  of  equation  (2.18).  This  equation  is  solved  first  with 
an  assumed  pressure  gradient  in  the  manner  usually  applied  to  the  external 
boundary-layer  problem,  for  brevity  referred  to  here  as  the  standard  problem; 
the  next  values  of  dp*  are  obtained  by  applying  Newton's  method.  For  example, 
rewritting  dp*  by  -6  and  equation  (2.28)  as 

*(BV)  =  f(x,nL)  -nL  (3.30) 

then  the  next  value  of  is  obtained  from 


Bv+1  .  ev - - 

3/3B[*(Bv)] 


v  =  0,  1 ,  2,  . . . 


(3.31) 


The  derivative  of 


with  respect  to  e  is  obtained  from  equation  (3.30),  i.e. 


fe  C4>(BV)]  =  fj[f(;,nL)]  (3.32) 

The  derivative  of  f  with  respect  to  6  is  obtained  by  solving  a  system  of 
variational  equations  to  be  described  later  in  this  subsection.  The  iteration 
process  is  repeated  until 

|Bv+1(x)  -  Bv+1(x)|  <  z}  (3.33) 

where  e-|  is  a  small  error  tolerance  which  may  be  specified. 

The  standard  problem  and  the  variational  equations,  described  in  detail 
below,  are  solved  with  the  two-point  finite-difference  method  developed  by 
Keller  [2],  The  application  of  the  standard  problem  to  two-  and  three- 
dimensional  boundary  layers  is  described  in  detail  in  reference  28  and  only 
a  brief  outline  is  provided  here. 
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In  accordance  with  Keller's  Box  scheme,  the  new  dependent  variables 
u(x,n)  and  v(x,n)  are  introduced  and  allow  equation  (2.18)  to  be  written  as 
the  first-order  system. 


f'  =  u 
u'  =  v 

(b  v).  +  l  fv  =  x(d£l+u^-v^ 

\dx  ax  ax/ 


(3.34a) 

(3.34b) 

(3.34c) 


On  0  <_  x  <_  Xj^j,  0  £  n  i  n„»  we  place  a  possibly  nonuniform  net: 


K  -  °» 

*n  =  Vi  +  kn 

n  =  1 ,  2,  ...»  N 

0 

(3.35) 

n0  -  0, 

"3  5  "3-1  *  h3 

4  =  1  >  2,  0  nj  = 

with  kn  and  h.  denoting  variable  distances  between  nodes  in  the  x  and  n 
directions.  The  finite-difference  forms  of  equations  (3.34a)  and  (3.34b),  with 
central -difference  quotients  and  averages  about  the  midpoint  (xn,  nj_i/2)  > 
then  can  be  written  as 


(fj 


)  * 

)  = 


Uj-l/2 

"J-l/2 


(3.36a) 

(3.36b) 


Similarly  equation  (3.34c)  Is  approximated  by  centering  about  the  midpoint 
(*n-l/2’  n j _ i / 2 )  a  rectangle  whose  mesh  widths  are  kn  and  hj: 


hj  C(b-,v)j  (b1v)”-1]  an(u2)J_1/2  +  J  +  an  ^fv^j-l/2 

+  an[fj-1/2vj-l/2  “  vj-l/2fj-l/2]  "  T j -1/2  “  2an6 


(3.36c) 
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where 


T 


n-1 

j-1/2 


n^(fv)j-l/2  ”  (u2)j-l/2 
.  -Vl/2 


d -  (biv)-1 

6  =  (p*)n_1  -  (p*)n 


] 


Equations  (3.36)  are  imposed  for  j  =  1,2,...,  J.  The  wall  boundary 
conditions  of  equation  (2.23b),  with  the  new  dependent  variables  defined  in 
equation  (3.34),  yield,  at  x  =  xn. 


*  0  uj  =  0  (3.37a) 

Similarly,  the  outer  boundary  conditions  of  equation  (2.24)  can  be  written 
as 

(u2)J  =  (u2)J_1  +  2b  (3.37b) 

If  (fj”\  uj »  vj"^)  an(1  6n  are  assumed  known  for  0  _<  j  <_  J ,  then 
equation  (3.36),  for  1  <_  j  £  J-l  and  the  boundary  conditions  of  equation 
(3.37),  yield  a  nonlinear  algebraic  system  of  3J+3  equations  in  as  many 
unknowns  (f!-,  Uj,  Vj).  The  system  can  be  easily  solved  by  block  elimination 
after  linearization  by  Newton's  method.  Details  of  Newton's  method  and  the 
block  elimination  method  for  the  momentum  equation,  including  the  corresponding 
computer  program  program,  are  provided  in  [28]. 

To  determine  the  value  of  3<j>/38,  in  equation  (3.31),  it  is  necessary  to 
know  3f/3B  at  (xn,n[_).  Thus,  equation  (3.36)  is  differentiated  with  respect 
to  B  and  leads  to  the  following  linear  difference  equations,  known  as  the 
variational  equations  for  equation  (3.36): 
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4. 


-  —  . -  ^-^-u*4**H* 


ii  iiftiia 


(3.38a) 


h>3  -g3-!>  =  UJ-l/2 

hJ1(Uj  ~  UJ-l'  =  Vj-l/2 

hJ,[(bl)jVJ  _  _“n(uJUj  +  uJ-lUj-l * 


-  (1  -  +  vjG5  -  -  vj-igj-1) 

+  ?«n[,J-l/2(6l+83-l>-tJ-WVi  +  V3-1)]- 


(3.38b) 

(3.38c) 

-2* 


Similarly,  the  boundary  conditions  of  equation  (3.37)  yield 


g"  =  0 


u"  =  0 
o 


u3  ■  i/u3 


(3.39) 

In  equations  (3.38)  and  (3.39),  the  following  notation  has  been  introduced: 


G  E  36  U  5  I?  V  s  f  (3.40) 

In  common  with  the  standard  problem,  equations  (3.38)  and  (3.39)  form  a 
tridiagonal  system,  with  3x3  blocks,  which  is  solved  by  the  block  elimination 
method  referred  to  above  and  described,  in  detail,  in  reference  28. 

3.3  Momentum  Equation  and  the  Mechul  Method  for  Internal  Flows 

In  contrast  to  the  nonlinear  eigenvalue  approach,  the  Mechul  method  solves 
directly  for  e  and  eliminates  the  need  for  iteration.  Thus  the  conservation 
equations,  in  transformed  variables,  again  contain  one  unknown  more  than  the 
number  of  equations.  Mass  continuity,  in  integral  form,  together  with  a 
knowledge  of  the  duct  geometry  provides  the  required  additional  information 
and  allows  explicit  solution  of  the  equation,  without  iteration,  and  with 
pressure  obtained  from  the  solution.  With  the  function  p*  treated  as  unknown 
and  with  the  help  of  the  y-momentum  equation  (-3p*/8y  =  0)  equation  (2.18) 
can  be  written 
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*'  s  u  (3.41a) 

«'  ■  v  (3.41b) 

(P*)'  =  0  (3.41c) 

(b1v) *  +  y  fV  =  x  +  u  v  (3.41d) 

\dx  ax  ax/ 


and  the  wall  boundary  conditions  (2.23b)  become: 


n=°  f  =  u  =  0  (3.42a) 

Since  as  y  -*■  6,  f  ^  n»  an  "edge"  boundary  condition  for  f  at  n  =  n 

oo 

can  be  written  with  equation  (2.28).  Then  the  boundary  conditions  at  the 
edge  become 

n  =  n-  f  s  ue^n»  “  +  u  =  ue  (3.42b) 

The  difference  equations  for  equation  (3.41a,b,c)  are  similar  to  those 
given  by  equations  (3.36a,b),  and  those  for  (3.41d)  are  similar  to  those  of 
equation  (3.36c),  except  that  B  is  now  variable;  that  is 


(3.43a) 

"SX  '“"-I’  ‘“3-1/2 

(3.43b) 

hi1(pj  -  pj-i' -  0 

(3.43c) 

^  +  {7*  an)*fv>j-l/2  _an[2|>J-)/2  + 

(u2)n 

1  j-1/2 

Vj-l/2fj-l/2  +  fj-l/2Vj-l/2]  =  Rj-l/2  (3*43d) 


*3-1/2  *  an[(fv)j-l/2  (u2)j-l/2  ~  ^j-l/2]  _hj1  [(blv)j_1  ~  (blv)j-l] 


2  (fv)j-l/2 
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Similarly  the  boundary  conditions  (3.42a,b)  become 


fj  =  u"  =  0, 
o  o 


“3  ■  “3* 


(3.44) 


The  solution  of  the  system  given  by  equations  (3.43)  and  (3.44)  is 
similar  to  that  for  the  standard  problem.  This  time  there  is  a  nonlinear 
algebraic  system  of  4J+4  equations  in  as  many  unknowns  (f!j,  uj,  v",  p!j). 

Again  the  system  can  be  solved  with  the  block  elimination  method  discussed 
earlier  after  linearization  by  Newton's  method. 

For  flows  with  negative  wall  shear,  it  is  necessary  to  make  approximations 
to  the  governing  equations  to  continue  the  calculations  past  the  separation 
point.  Here  we  use  the  approximation  first  suggested  by  Reyhner  and  Flugge-Lotz 
[32J.  This  approximation  referred  to  as  FLARE  by  Williams  [24],  neglects  the 
u(3u/3x)  term  in  the  region  of  negative  u-velocity,  u  <  0.  All  inverse 
boundary-layer  procedures,  including  the  Mechul  method  described  above,  use 
this  approximation  for  regions  of  separated  flow.  Once  a  solution  has  been 
obtained  with  this  approximation,  however,  the  assumption  can  be  relaxed  to 
allow  solutions  of  the  complete  boundary-layer  equations.  For  details,  the 
reader  is  referred  to  [11]. 


3.4  Energy  Equation  for  Internal  Flows 

The  method  used  to  obtain  solutions  to  equation  (2.19)  is  similar  to  that 
described  in  section  3.2  except  that  the  iteration  technique,  necessitated  by 
the  presence  of  pressure  gradient,  is  no  longer  required.  Again,  a  new 


dependent  variable  w(X,n)  is  defined,  i.e. 


g 1  =  w 


(3.45a) 


and  equation  (2.19)  rewritten  as: 


(b2w) 1  +  1  Pr  fw  =  Pr  x  /u  -  w  — \  (3.45b) 

\  3x  3 x  / 
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Equation  (3.45a)  is  centered  at  (J?n,nj_1/2)  and  equation  (3.45b)  at 
(*n-l/2’  nj-l/2^  to  obtain  the  finite-difference  equations. 

9j  "  9J-1  "  /  (wj  +  wj-l  >  =  0  0.46a) 

<51 )jwj  *  (s2>jwj-l  *  s3<9j  +  9j-1>  '  <r1>j  0.46b) 

Here 


ki*j =  — * Pr  [i  fj- 


+  „  /fn  *n-l  / 

1/2  °n  j-1/2  j-1/2  _ 


,  ,  (K2>?  ,  (b,)" 

U2>J - (C1  >J 

k3>j  -  -2pr  “„ujn:l/2 


(3.47) 


(r)>3  ■ 


-1 

1/2 


+  j-1/2  f j - 1 / 2  ^ w j - 1 / 2] 


The  solution  of  the  linear  system,  subject  to  the  general  wall  boundary 
conditions 

ao9o  +  alw2  =  Yo  (3.48a) 

and  to  the  general  "edge"  boundary  conditions 

So9J  +  elwJ  '  Y1  (3.48b) 

can  be  obtained  by  the  same  block  elimination  method  used  with  the  momentum 
equation. 
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r 


n 


The  system  of  equations  (3.46)  and  (3.47)  has  a  block  tridiagonal 


structure  which 

is  shown 

in  matrix-vector  form 

r  go 

wo 

9J 

w.  g 

b.c.  eq. (3.48a) 

1 - 

i  a 
l  0 

“1 

i  0 

1 

"i"! 

1 

eq.  (3.46a) 

1 

1 

I  -1 

-hl 

2 

1 

1 

1  1 

1 

a  i 

2 

eq.  (3.46b) 

r<c3>j 

i 

<*;r 

1 

l  0 

i _ 

0 

i 

i  -i 
j _ 

_hj+l  i 

2  1 
i 

!  (t3>J 

(«2)jJ  (4 

b.c.  eq. (3.48b) 

1 

1  0 

l 

0  I  3 

J _ 

— i — 

w, 


-h. 


+1 


«rl>o\ 


‘(r2>o / 

'(rlV 


'(r2,jV 

'hV 

(r2)J, 


(3.49) 


This  equation  can  also  be  written  as  in  (3.25)  and  may  be  solved  by  the  pro¬ 
cedure  described  following  equation  (3.25). 

3. 5  Comments  on  the  Solution  Algorithm 

The  Box  method  allows  nonuniform  net  spacings  in  the  streamwise  direction 
and  across  the  boundary  layer.  The  use  of  nonuniform  net  spacing  across  the 
layer  is  essential  for  turbulent  flow  calculations.  This  can  be  achieved  by 
usinq  a  grid  which  has  the  property  that  the  ratio  of  lengths  of  any  two 
adjacent  intervals  is  a  constant;  that  is,  hj  =  Kh j _i .  The  distance  to  the 
j-th  line  is  given  by  the  following  formula: 

nj  =  hl(KJ  ”  D/(K  “  T)  K  >  1  (3.50) 
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There  are  two  parameters:  h-j ,  the  length  of  the  first  An  step,  and  K, 
the  ratio  of  two  successive  steps.  The  total  number  of  points  J  can  be 


calculated  from  the  expression, 


ln[l  +  (K  -  1  )(nj#/h1 )] 


(3.51) 


For  further  details,  see  Cebeci  and  Bradshaw  [28] 
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4.0  RESULTS 


The  present  results  section  has  been  prepared  in  six  subsections.  The 
first  considers  laminar  internal  flows  and,  in  view  of  the  greater  importance 
of  turbulent  flow  and  the  considerable  attention  previously  devoted  to  this 
topic,  is  not  extensive.  It  serves,  however,  to  provide  confidence  in  the 
abilities  of  the  numerical  method  and  to  demonstrate  its  efficiency.  The 
eddy-viscosity  and  turbulent  Prandtl -number  assumptions  are  described  in 
subsection  4.2  and  briefly  compared  with  available  alternatives.  Turbulent 
internal -flow  results,  obtained  with  the  assumptions  of  subsection  4.2  and 
with  the  numerical  method  of  subsections  3.2  and  3.4  are  presented  in  sub¬ 
section  4.3  and  compared  with  previous  calculations  and  measurements.  New 
internal-flow  results  are  presented  in  subsection  4.4  to  quantify  the 
influence  of  Reynolds  number  on  the  Stanton  number  distribution  for  pipe  flow, 
including  regions  of  laminar,  transition  and  turbulent  flow.  Separated  flows 
are  considered  in  subsections  4.5  and  4.6  where  the  Mechul  approach  of  subsection 
3.3  is  applied  to  laminar  and  turbulent  external  flows,  respectively.  The 
internal -flow  results  are  based  on  those  of  reference  8  and  the  separated-flow 
results  of  reference  9. 

4.1  Laminar  Internal  Flow 

It  is  difficult  to  find  useful  practical  applications  for  results  relating 
to  the  entrance  region  of  a  pipe  flow,  where  the  initial  profile  is  presumed 
uniform.  Nevertheless,  a  large  number  of  computational  investigations  have 
been  reported  and  a  sample  are  referred  to  in  table  1.  Corresponding  experi¬ 
ments  have  also  been  reported  largely  because  of  the  interest  stemming  from 
transpiration  viscometers  of  finite  inlet  length.  In  practice,  the  initial 


profile  depends  upon  the  shape  of  the  entrance  to  the  tube  and  on  the 
expansion  ratio,  as  is  witnessed  by  the  experimental  results  referred  to  in 
table  1.  It  may  be  assumed  that  for  a  sharp-edged  entrance  normal  pressure 
gradients  and  small  regions  of  flow  separation  will  exist  and  may  invalidate 
the  boundary-layer  assumptions. 

Nevertheless,  and  to  provide  a  sample  laminar-duct-flow  calculation  which 
may  be  compared  with  previous  results,  the  laminar  entrance-flow  problem  has 
been  considered  and  solutions  obtained.  A  sample  of  the  results  are  reproduced 
by  Cebeci  and  Bradshaw  [28]  and  may  be  supplemented  by  figure  6  which  presents 
the  calculated  distribution  of  nondimensional  pressure  drop  from  the  tube 
entrance,  where  uniform  velocity  has  been  assumed,  to  the  downstream  parabolic 
profile  region.  The  solutions  of  Langhaar  [33]  and  Campbell  and  Slattery  [34] 
are  also  shown  together  with  the  Poiseuille  flow  asymptote.  Experimental  data 
have  also  been  reported,  for  example  by  Pfenninger  [35],  Shapiro,  Siegel  and 
Kline  [36]  and  by  Whitelaw  [37]  and,  as  can  be  seen,  for  example,  from  the 


Table  1. 


Kinetic  energy  correction  factor  h  defined  by 


P  ~  P0 

- - 2*-  =  64x**  +  2h 


Calculated 


Present  Investigation  1.123 

Langhaar  [33]  1  J4 

Campbell  and  Slattery  [34]  1*09 

Measured 

Schiller  [38]  1.058  -  1.225 

Rieman  [39]  1.110-1.134 

Swindells  et  al  [40]  1.12  -  1  17 
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Figure  6.  Pressure  drop  for  laminar  flow  in  a  circular  pipe. 


paper  of  Campbell  and  Slattery,  are  subject  to  scatter  but  tend  to  support  the 
present  results.  Since  viscometry  is  a  likely  application  of  results  to  the 
entrance-length  problem,  table  1  presents  calculated  and  measured  values  of  the 
so-called  "kinetic-energy  correction  factor." 

An  example  of  heat-transfer  calculations,  for  flow  in  a  two-dimensional 
plane  channel  with  developing  momentum  and  thermal  boundary  layers  with 
constant  heat  flux  are  presented  in  figure  7.  The  results  may  be  compared 
with  the  previous  calculations  of  Siegel  and  Sparrow  [41]  and  Naito  [42], 
which  have  been  represented  by  a  number  of  discrete  points  taken  from  the 
plotted  results,  and  are  in  excellent  agreement  except  at  the  lowest  values  of 
the  Prandtl  number.  The  more  recent  results  of  Bhatti  and  Savery  [43],  which 
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also  made  use  of  the  Karman-Pohl hausen  method  to  solve  the  energy  equation,  are 
not  displayed  but  are  in  general  agreement  although  somewhat  lower  than  present 
and  previous  results;  the  difference  probably  stems  from  the  velocity-profile 
assumption  of  Bhatti  and  Savery. 

In  the  above  calculations  of  laminar  flow,  the  standard  calculations  per¬ 
formed  upstream  of  the  coalescence  were  obtained  by  using  approximately  10  to 
15  x-stations.  In  the  region  downstream  of  boundary-layer  coalescence,  we  have 
used  approximately  15  to  20  x-stations.  These  large  forward  steps,  which 
contributed  greatly  to  the  efficiency  of  the  present  calculation  method,  were 
made  possible  by  the  efficiency  and  stability  of  the  finite-difference  scheme 
and,  in  the  upstream  region,  by  the  use  of  transformed  coordinates. 


4.2  Physical  Assumptions  for  Turbulent  Flow 

For  turbulent  flows  in  ducts,  there  are  essentially  three  distinct  regions. 
Near  the  entrance,  the  shear  layers  develop  in  the  region  close  to  the  solid 


Figure  7.  Local  Nusselt  number  distributions  for  laminar  flow  in  a  plane  channel 
with  constant  wall  heat  flux. 


wall  and  the  core  velocity  increases  continuously  due  to  the  growth  in 
displacement  thickness.  The  flow  in  this  region,  which  may  be  called  the 
displacement-interaction  region  [43],  is  physically  no  more  than  the  boundary- 
layer  flow  in  a  mild  favorable  pressure  gradient. 

As  the  flow  moves  downstream,  the  shear  layers  grow  and  finally  meet  on 
the  centerline.  The  shear-layer  interaction  increases  in  strength  until  a 
steady  state  of  mutual  eddy  intrusion  and  fine-scale  mixing  is  achieved,  where 
the  flow  is  fully  developed. 

The  eddy-viscosity  concept  is  used  here  to  model  the  Reynolds  shear  stress, 
-pu'v1  and,  in  the  displacement  interaction  region,  the  formulation  is  that  of 
reference  18: 


Inner  layer 


ej  =  T<y[l  -  exp(-y/A)]}‘ 


e .  <  e 

1  -  0 


(4.1a) 


Outer  layer 


eQ  =  0.0168 


Here 


k  =  0.40, 


OD 

/  (ue  -  u)dy 


A  =  26v/ut, 


(4.1b) 


UT  =  y/‘rv/P  (4-2) 


In  the  fully  developed  region,  we  use  the  well-known  mixing-length 
formula  of  Nikuradse  and  express  the  eddy  viscosity  for  the  entire  layer  by 


au 

ay 


e  *  j*Cl  -  exp(-y/A)]j2 


(4.3) 


where 

i  =  rQ[0.14  -  0.08(1  -  r/rQ)2  -  0.06(1  -  r/ro)4]  (4.4) 

To  use  equation  (4.4)  for  two-dimensional  duct  flows,  we  simply  set  rQ  =  h 
to  get 

1  -£-•  y/h 
0 

For  the  region  between  the  displacement  interaction  and  fully  developed 
regions,  the  eddy-viscosity  formula  is  given  by  the  following  expression  that 
combines  equation  (4.1)  and  (4.4) 


e 


1  -  exp 

(x-x)l 

U 

xr 

> 

L  o  J 

(4.5) 


with  ej  denoting  the  formula  given  by  equation  (4.1),  ep  the  formula  by 
equation  (4.3),  and  with  x  denoting  an  empirical  constant. 

In  cases  where  the  location  of  the  onset  of  transition  is  known,  the 
flow  may  be  determined  through  the  transition  region  and  into  the  fully 
turbulent  region  by  multiplying  the  expressions  for  e  by  an  intermittency 
factor  Ytr,  defined  by 


Ytr  = 


1  -  exp 


1  R0 

T20IT  K" 


66 

tr 


(4.6) 


For  calculations  requiring  solutions  of  the  energy  equation  for  turbulent 
flows,  a  turbulent  Prandtl  number  may  be  specified.  Here  we  use  the  formula 
described  in  [1].  According  to  that  formula,  the  turbulent  Prandtl  number 
Prt,  is  computed  from 


Pr 


> . 40[ 1  - 
l.44[l  - 


-  exp 


exp 


(4.7) 
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with 


.  j 

B  =  B+v/u_,  B+  =  - 4™-  V  Ci  (log  Pr) 

(Pr)1^  Y 


i-1 


(4.8) 


C]  =  34.96,  C2  =  28.79,  C3  -  33.95,  C4  =  6.33  and  Cg-  -1.186 


4.3  Turbulent  Internal  Flow 

For  the  convenience  of  the  reader,  a  sample  of  the  results  obtained  for 
turbulent  flow  are  presented  in  the  following  paragraphs.  The  capabilities  of 
the  turbulence  model  and  its  relationship  to  alternatives  are  discussed  after 
the  results  and  the  comparisons  with  measurements  have  been  introduced. 

Most  of  the  turbulent  flows  considered  here  relate  to  plane  channels,  but 
to  preserve  symmetry  with  section  4.1,  it  is  useful  to  start  with  pipe-flow 
results.  The  measurements  of  Barbin  and  Jones  [45]  provide  an  example  of  a 
developing  pipe  flow  which,  in  the  near  fully-developed  condition,  accords 
with  the  classical  measurements  of  Laufer  [46].  It  should  be  emphasized  that 
the  "fully-developed"  conditions  are  achieved  asymptotically  and,  at  their 
diameter  Reynolds  number  of  3.88  x  10^,  Barbin  and  Jones  concluded  that  it 
had  not  been  achieved  at  40.5  diameters.  The  pressure-drop  results  of  figure  8, 
which  are  perhaps  the  least  sensitive  check,  would  suggest  that  the  flow  was 
very  close  to  the  asymptote.  The  computed  nondimensional  centerline  velocity 
distribution  is  also  shown  on  figure  8a  and,  like  the  pressure-drop  data,  is 
in  close  accord  with  measurements  although  in  the  region  where  the  boundary 
layers  fill  the  pipe,  they  show  a  more  rapid  change  in  shape  than  the 
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measurements.  This  slight  discrepancy  between  calculations  and  measurements 
can  be  attributed  to  the  form  of  the  eddy-viscosity  hypothesis  and  is  not 
reflected  in  the  profiles  of  figure  8b.  The  calculations  began  with  the 
measured  profile  at  x/D  of  1.5. 

The  turbulent  flow  in  two-dimensional  channels  has  been  investigated  by 
many  authors  and,  together  with  the  pipe-flow  investigations,  has  added  a 
great  deal  to  an  understanding  of  turbulent  flows  in  general  and  near-wall 
flows  in  particular.  The  contributions  of  Laufer  [47]  and  Comte-Bellot  [48] 
are  perhaps  best  known  although  the  more  recent  paper  by  Hussain  and  Reynolds 
[49]  uses  the  largest  aspect  ratio  of  18:1;  this  is  not  to  say  that  the 
resulting  flow  is  necessarily  the  most  two-dimensional.  In  addition,  isothermal 
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Figure  8.  Comparison  of  present  calculations  with  Barbin  and  Jones'  pipe-flow 
measurements,  (a)  Centerline  velocity  and  pressure  drop,  (b)  Veloc¬ 
ity  profiles. 


two-dimensional  plane  channel  results  have  been  reported,  for  example  by 
Clark  [50],  Byrne,  Hatton  and  Marriott  [51]  and  Dean  [52],  The  sample  of 
calculated  results  shown  on  figures  9  to  11  has  been  presented  to  display 
a  range  of  properties  as  well  as  to  allow  comparison  with  the  results  of 
three  experiments. 

Figure  9  compares  calculated  velocity  profiles  with  the  measurements  of 
Comte-Bellot.  The  calculations  began  with  the  measured  profile  at  x/D  of 
20  and,  as  can  be  seen,  are  in  very  close  accord  with  the  measurements. 
Distributions  of  displacement  and  momentum  thickness,  corresponding  to 
the  experiment  of  Byrne  et  al .  [51]  are  presented  in  figure  10,  and  it  is 


Figure  9.  Comparison  of  calculated  velocity  orofiles  with  the 
measurements  of  Comte-Bellot. 
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Figure  10.  Computed  displacement  and  momentum  thickness  distributions  in  a  plane 
channel;  comparison  with  the  measurements  of  Byrne  et  al . 

clear  that  the  calculated  values  are  within  experimental  scatter.  Dean's  data 
[50]  for  velocity  profile  centerline  velocity  and  wall-shear  stress  are  pre¬ 
sented  in  figure  11;  the  centerline  velocities  are  very  well  represented  and 
the  wall-shear  stress  values  are  within  7%  of  measurements. 

The  results  of  figures  9  to  11,  in  general,  show  that  the  present  cal¬ 
culation  method  is  able  to  represent  turbulent  internal  flows  with  adequate 
accuracy  for  engineering  purposes.  Small  discrepancies  do,  however,  exist  and 
are  partly  attributable  to  the  present  eddy-viscosity  assumption  and  partly 
due  to  measurement  imprecision  and  deviations  from  two  dimensionality.  In  the 
flow  upstream  of  the  coalescence  of  the  boundary  layers  on  the  two  sides  of  the 
channel,  the  calculations  and  experiments  are  in  excellent  agreement  apart, 
perhaps,  from  the  near-wall  region  of  some  of  the  velocity  profiles  of  Comte- 
Bellot.  In  this  region,  the  inner-layer  expression  for  eddy  viscosity 
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Comparison  of  present  calculations  with  Dean's  plane  channel-flow 
measurements,  (a)  Centerline  velocity  and  shear  stress,  (b) 
Velocity  profiles,  (c)  Centerline  velocity  with  different  values 
of  the  parameter  X. 


dominates  the  calculated  results  which  are  strongly  dependent  on  th'.  value  of 
the  von  Karman  constant,  k.  The  value  of  k  used  in  this  region  was  0.4.  As 


pointed  out  by  Hussain  and  Reynolds,  however,  this  value  with  an  A  of  26 
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leads  to  lower  values  of  em>  in  the  near-wall  region,  than  their  experiments. 
It  appears  that  it  also  results  in  lower  values  than  the  experiments  of  Comte- 
Bellot  although  it  accords  with  the  results  of  Laufer. 

The  calculated  results  for  the  data  of  Barbin  and  Jones,  for  pipe  flow, 
also  suggest  that  a  higher  value  of  k  could  be  used:  the  profiles  of  figure 
8b  suggest  this  and  it  is  brought  out  even  more  clearly  by  the  selection  of 
related  profiles  presented  by  Cebeci  and  Bradshaw.  The  centerline  velocity 
distributions  presented  on  figures  8a,  11a  and  c,  for  a  pipe  and  channel, 
respectively,  indicate  excellent  predictions  up  to  the  region  of  coalescence. 
In,  and  downstream  of  this  region,  the  agreement  with  Dean's  results  is 
excellent  for  a  value  of  A  of  20  and  less  so  with  the  results  of  Barbin  and 
Jones.  It  is  clear  from  figure  11c,  however,  that  a  finite  value  of  A  is 
necessary  and  computational  experiments  have  suggested  the  value  of  20  used 
for  present  calculations. 

The  values  of  wall  shear  stress,  indicated  in  figure  11a,  are  consistent 
with  the  values  of  k  and  A+  and  with  the  prediction  of  lower  velocity 
values  in  the  near-wall  region.  Adjustment  of  these  constants  would  allow  the 
values  of  shear  stress  to  accord  with  measurements  without  altering  the  results 
of  figure  lib.  The  results  obtained  with  the  interacting  version  of  Bradshaw's 
stress  model  (see  Bradshaw,  Dean  and  McEligot  [16])  are  less  satisfactory  for 
the  centerline  velocity  distribution  but  more  so  for  wall  shear  stress,  partly 
due  to  the  use  of  0,41  as  the  von  Karman  constant. 

Both  Na  and  Habib  [13]  and  Stephenson  [53]  have  obtained  calculated  distri 
butions  of  Nusselt  number  as  a  function  of  Reynolds  number  for  fully-developed 
turbulent  flows  In  round  pipes  and  for  a  Prandtl  number.  As  indicated  by 
Habib  and  Na  [14],  the  results  are  essentially  independent  of  the  temperature 
boundary  condition  and  this  is  supported  by  the  present  results  which  are  shown 
on  figure  12  for  a  laminar  Prandtl  number  of  0.72.  As  expected,  the  resulting 
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Figure  12.  Computed  local  Nusselt  number  for  fully-developed  pipe  flow  with 
constant  wall  heat  flux  and  constant  wall  temperature;  comparison 
with  experimental  data.  For  references  see  [13]. 

curve  is  very  similar  to  those  of  Na  and  Habib.  For  a  bulk  Reynolds  number  of 
3.88  x  10®.  the  value  of  Stanton  number  is  1.85  x  TO"3  which  may  be  compared 
with  Stephanson's  value  of  2.0  x  10~3.  (The  corresponding  results  for  fully- 
developed,  two-dimensional,  plane  channel  flow  with  a  bulk  Reynolds  number  of 
1.92  x  105  are  2.19  x  10"3  and  2.47  x  10"3.) 

The  use  of  the  present  eddy-viscosity  and  turbulent  Prandtl  number  formu¬ 
lations  is  clearly  successful  in  that  they  allow  the  flows  examined  in  this 
section  to  be  very  well  represented.  It  is  possible  to  make  use  of  simple  or 
more  complicated  approaches.  The  present  eddy  viscosity  is  probably  the  sim¬ 
plest  which  can  be  used  to  represent  external  and  internal  boundary  layers 
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and  still  allow  the  possibility  of  calculations  from  transition,  as  demonstrated 
in  the  next  section.  It  cannot,  however,  be  expected  to  represent  flows,  such 
as  large  regions  of  separation,  where  the  flow  is  far  from  equilibrium  and 
transport  of  turbulence  properties  is  important.  On  the  other  hand,  the  use 
of  turbulence  models  which  allow  for  transport,  for  example,  of  turbulence 
energy  are  unnecessary  for  the  present  boundary- layer  flows. 


4.4  Internal  Flow  with  Transition 

Figure  13  is  concerned  with  a  flow,  reported  by  Hall  and  Khan  [54],  which 
developed  in  a  round  pipe  from  a  laminar  state,  through  transition  to  fully- 
developed  turbulent  flow.  Measurements  were  obtained  with  uniform  wall  heat 
flux  and  temperature  applied  from  close  to  the  entrance.  The  figure  allows 
comparison  between  measurements  and  calculations  at  pipe  Reynolds  numbers  of 
31,100  and  33,000  and  for  the  two  thermal -boundary  conditions.  The  flow  was 
assumed  to  undergo  transition  at  x  =  4D  where  the  Stanton  number  is  minimum 
as  shown  in  the  experimental  data.  As  can  be  seen,  the  agreement  is  good  and 
the  influence  of  the  thermal  boundary  condition  much  greater  in  the  region  of 
transition  than  in  the  laminar  and  turbulent  regions.  The  results  of  Hall  and 


(a) 


Figure  13.  Computed  local  Stanton  number  in  developing  pipe  flow;  comparison 
with  measurements  of  Hall  and  Khan  for:  (a)  Constant  heat  flux, 
(b)  Constant  wall  temperature. 
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Khan  show  that  the  influence  of  the  thermal  boundary  condition  diminishes  with 
increasing  Reynolds  number  and  present  calculations  support  this  finding.  Thus, 
the  uniform  wall  heat-flux  results  of  figure  14  accord  with  uniform  wall- 
temperature  results  for  values  of  RD  greater  than  4  x  104  to  better  than 
5"  of  the  local  value.  The  results  of  figure  14  were  obtained  with  the 
assumption  that  the  transition  occurred  at  a  value  of  momentum-thickness 
Reynolds  number  of  320. 

4.5  Laminar  Separated  Flows 

For  laminar  flows  with  separation,  four  separate  external  flows,  which 
allow  the  displacement  thickness  distribution  to  be  specified,  have  been 
considered.  These  include  the  laminar  separated  flows  of  Williams  [25], 


io4 

2  x  104 
4  x  104 
10S 

2  x  10S 
4  x  105 


Figure  14.  Local  Stanton  number  distributions  in  developing  pipe  flow  as  a 
function  of  Reynolds  number  (Pr  =  0.72  and  Retr  =  320). 
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Carter  [22]  and  Briley  [26].  In  the  last  case,  the  freestream  boundary 
condition  corresponded  to  a  linearly  decreasing  external  velocity  distribution 
followed  by  a  constant  velocity  and  resulted  In  separation  and  reattachment. 
Figure  15  allows  a  comparison  of  values  of  local  skin-friction  coefficient 
Cf  calculated  with  the  present  method  with  those  obtained  by  Briley  from 
the  steady,  two-dimensional  form  of  the  Navler-Stokes  equations.  Here  cf  is 
defined  by 


cf 


Tw  m  2f"(0) 

1/2  pu^  ^ 


(4.9) 


and  c^/Rl  represents  2f"(0).  These  calculations  were  made  from  the  dis¬ 
placement  thickness  distribution  (see  figure  15)  deduced  from  the  Navler-Stokes 
solutions.  In  general,  the  agreement  is  very  good  and  the  small  discrepancy 
may  be  associated  with  the  difficulty  of  reading  the  input  5*(x)  distribution 
from  the  graph  presented  by  Briley.  The  Navler-Stokes  solutions  were  obtained 
for  R|_  *  2.08  x  10^  and  required  45  minutes  on  a  UNIVAC  1108.  The  results 
obtained  by  the  Mechul  function  approach  of  subsection  3.3  required  approximately 
10  seconds  on  a  CDC-6600  computer. 


Figure  15.  Calculated  local  skin-friction  coefficient  distribution  for 
separating  and  reattaching  flow  computed  by  Briley  [26], 
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We  next  consider  the  two  laminar  flows  with  separation  and  reattachment 
computed  by  Carter  [22].  These  flows  have  displacement-thickness  distributions 
6*(x)  given  by 


1.7208/x 

6*(x)  =  l  a1  +  a,(x  -  x,)  +  a,(x  -  x,)2  +  a^(x  -  x,) 


with 


'1 


1'  T  °4' 
3 


'1 


a,  +  a,(x  -  x,)  +  a.  (x  -  x,) 


a1  =  1.7208 


1.0  <  x  <  Xj 

XT  1  x  1  x2 


x2  —  x  -  x3 


a2  =  (0.5)0.7208)/^ 

a3  =  (0.5/a^  («*ax  —  a-j )  -  4a2] 

a4  =  z/il  £i,/2a2  -  (s^, -a,)! 


a,  =  5* 

1  max 


V-’/Wr'Si'Vj'1 

s4  ’  ,/42  cj(4tox-53>  +  VS'3 


—  x2  x^  *  &2  —  x^  x2  ,  X-|  —  1 .065 .  x2  =  1 .35  j 


1.884 


(4.10) 


The  first  flow,  referred  to  as  Case  A,  has  6*  =  5.6,  and  the  second  flow 

max 

referred  to  as  Case  B,  has  6*ax  =  8.6.  The  two  displacement  thickness 
distributions  are  shown  in  figures  16  and  17. 

Comparison  of  present  results  with  those  of  Carter  [22]  is  shown  in 
figures  16,  17,  and  18.  The  present  calculations  were  started  at  x  =  0  by 
solving  the  governing  equations  in  transformed  variables  for  the  standard 
problem  and  then  at  x  =  1  the  Mechul -function  method  of  subsection  3.3  was  used 
to  solve  the  inverse  problem  with  the  equations  expressed  in  physical  variables. 
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Figure  17.  Local  skin-friction  distribution  for  Case  B 


Figure  18.  Streamline  pattern  in  separation  bubble  for  Case  B. 

As  can  be  seen,  the  present  results  agree  well  with  those  of  Carter  which 
made  use  of  the  Reyhner  and  Flugge-Lotz  approximations.  Our  skin-friction 
results  are  in  better  agreement  with  Carter's  results.  His  Reyhner  and 
Flugge-Lotz  approximation  was  slightly  different  from  that  considered  here 
and  may  account  for  the  disagreement  with  the  present  results.  The  computer 
time  associated  with  the  results  was  again  of  the  order  of  10s  on  a  CDC  6600. 
At  almost  all  x-stations  the  convergence  rate  was  quadratic  and  required  only 
two  or  three  iterations,  including  regions  of  separated  flow. 

Figure  19  shows  the  results  for  another  laminar  flow,  with  prescribed 
displacement  thickness  distribution,  previously  computed  by  Williams  [25]. 

The  6*-distribution  is  given  by 

6*  =  0.6479{1  +  9  exp[-0.0625(x  -  14)2]} 
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(a)  Wall  shear  distribution. 


(b)  External  velocity  distribution. 
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(c)  Velocity  profiles  near  and  In  the 
reverse  flow  region 

Figure  19.  Comparison  of  calculated  results  for  the  separated  flow  for  which 
displacement- thickness  distribution  Is  given  by  equation  (4.10). 


In  this  case  the  flow  starts  at  the  stagnation-point  and  separates  at  approxi¬ 
mately  x  -  10.  Both  calculations  were  made  by  the  Reyhner  and  Flugge-Lotz 
approximation.  As  in  the  previous  cases,  the  solutions  converged  quadratically 
requiring  only  two  or  three  iterations,  including  regions  of  separated  flow. 

4.6  Turbulent  Separated  Flows 

Two  separate  turbulent  flows  have  been  considered.  These  flows  did  not 
have  flow  separation  and  were  computed  with  the  eigenvalue  and  the  Mechul- 
function  methods  to  check  the  predictions  of  both  methods  for  attaching  flows. 
The  two  sets  of  results  were  identical.  The  application  of  the  method  for 
turbulent  flows  with  separation  is  currently  under  examination. 

The  two  attached  turbulent  flows  were  those  measured  by  Schubauer  and 
Spangenberg  (see  reference  27)  and  labeled  as  4400  and  4800  at  the  1968 
Stanford  Conference.  Flow  4400  has  a  strong  adverse  pressure  gradient  and 
flow  4800  has  a  mild  adverse  pressure  gradient.  Figures  20  and  21  present 
the  results  for  these  two  flows.  The  calculations  were  made  first  by  using 
the  standard  procedure:  in  this  case  the  external  velocity  distribution  was 
specified  and  the  boundary- layer  parameters  were  computed.  Next,  the  calcula¬ 
tions  were  made  with  both  of  the  inverse  procedures:  the  measured  displacement 
thickness  distribution  was  specified  and  the  external  velocity  distribution 
computed  along  with  other  boundary-layer  parameters.  As  the  results  show, 
the  agreement  with  experiment  is  very  good;  the  difference  between  calculated 
and  experimental  velocity  distributions  is  small,  and  the  boundary-layer 
parameters  computed  by  standard  and  inverse  procedures  agree  with  each  other. 


Figure  20.  Results  for  flow  4400.  Figure  21.  Results  for  flow  4800. 

(Inverse  results  were  (Inverse  results  were 

obtained  by  using  both  the  obtained  by  using  both  the 

nonlinear  eigenvalue  method  nonlinear  eigenvalue  method 

and  the  Mechul -function  and  the  Mechul -function 

method  and  are  identical.)  method  and  are  identical.) 


5.0  CONCLUDING  REMARKS 


The  results  of  the  previous  sections  allow  an  evaluation  of  the 
calculation  procedure  in  terms  of  the  finite-difference  method  and  the 
physical  assumptions.  Relevant  comments  have  been  included  in  section  4  and 
are  summarized  here. 

Comparison  of  computer  run  times  for  internal-flow  calculations  with  those 
of  alternative  numerical  methods  shows  that  the  present  formulation  is  advan¬ 
tageous  for  the  boundary- layer  problems  which  it  is  intended  to  solve.  The 
reduction  of  the  second-order  equations  to  first-order,  the  introduction  of 
transformed  coordinates,  and  the  use  of  the  "Box"  method  of  finite-differencing 
ensures  that  the  efficiency  is  coupled  with  second-order  accuracy.  A  typical 
computation  time  for  a  case  using  41  nodes  across  the  duct,  and  35  axial  sta¬ 
tions  down  the  duct  was  less  than  2  seconds  on  an  IBM  370/175.  The  chosen 
turbulence  model  is  near  optimum  for  the  range  of  flow  problems  considered  and 
has  been  shown  to  perform  well,  see  for  example,  Cebeci  and  Smith  [1]  and 
Cebeci  and  Bradshaw  [28],  for  a  much  wider  range  of  boundary-layer  problems 
including  two-dimensional  and  three-dimensional  external  flows. 

The  separated  flow  results  show  that  the  present  numerical  scheme,  with 
the  Reyhner  and  Fliigge-Lotz  approximation  [32]  can  be  used  to  represent  small 
regions  of  separated  flow.  The  use  of  boundary-layer  equations  to  represent 
separated  flows  presumes  that  the  normal  pressure  gradient  and  longitudinal 
diffusion  are  small  and  the  range  of  successful  applicability  remains  to  be 
determined.  Nevertheless,  it  has  been  shown  that  the  more  efficient  and  accu¬ 
rate  solution  of  boundary-layer  equations  will  compensate  for  the  approximations 
over  a  range  of  flows  with  small  regions  of  recirculation.  The  more  expensive 
solution  of  the  elliptic  equations,  including  the  normal  pressure-gradient  and 
longitudinal-diffusion  effects,  becomes  Increasingly  so  as  the  number  of  grid 


points  Is  Increased  to  make  a  finite-difference  solution  converge  to  the  solu¬ 
tion  of  the  equations,  see  for  example  reference  10. 

The  following  more  important  conclusions  may  be  extracted  from  the  present 
study: 

1.  The  present  numerical  scheme  allows  the  solution  of  boundary-layer 
equations,  for  a  wide  range  of  boundary  conditions,  with  efficiency 
and  accuracy.  There  are  no  significant  difficulties  in  applying  it 
to  internal -flow  problems. 

2.  The  present  eddy-viscosity  and  turbulent  Prandtl  number  formulations, 
previously  used  in  a  wide  range  of  external  boundary- layer  flows,  has 
been  shown  to  be  equally  satisfactory  for  a  range  of  internal-flow 
problems. 

3.  The  Internal-flow  problems  considered  include  incompressible,  two- 
dimensional  laminar  and  turbulent  pipe  and  plane-channel  without  and 
with  heat  transfer.  Provided  the  location  of  the  onset  of  transition 
Is  known,  flow  problems  of  this  type  including  compressible  flows  can 
readily  be  solved  with  the  present  method. 

4.  The  nonlinear  eigenvalue  method  used  to  solve  the  internal-flow  prob¬ 
lems  can  be  replaced  by  the  Mechul -function  method  which  can  also  be 
used  to  solve  separated-flow  problems.  The  present  calculation 
demonstrate  that  the  Mechul -function  method,  with  a  Reyhner  and 
Fliigge-Lotz  approximation,  conveniently  leads  to  results  which  are 

in  close  agreement  with  previous  calculations  based  on  solution  to  the 
Navler-Stokes  equations.  Since  the  present  procedure  is  considerably 
less  expensive  than  solutions  of  the  Navler-Stokes  equations,  the 
range  of  applicability  should  be  established. 
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